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ABSTRACT 

We have derived the star formation history (SFH) of the blue compact dwarf galaxy I Zw 18 through 
comparison of deep HST/ACS data with synthetic color magnitude diagrams. A statistical analysis 
was implemented for the identification of the best-fit SFH and relative uncertainties. We confirm 
that I Zw 18 is not a truly young galaxy, having started forming stars earlier than ~1 Gyr ago, 
and possibly at epochs as old as a Hubble time. In I Zw 18' s main body we infer a lower limit of 
w 2 x 1O 6 M0 for the mass locked-up in old stars. I Zw 18' s main body has been forming stars very 
actively during the last ~10 Myr, with an average star formation rate (SFR) as high as » IMq/ut 
(or k 2 x 10~ 5 MQyr~ l per 2 ) . On the other hand, the secondary body was much less active at these 
epochs, in agreement with the absence of significant nebular emission. The high current SFR can 
explain the very blue colors and the high ionized gas content in I Zw 18, resembling primeval galaxies 
in the early Universe. Detailed chemical evolution models are required to quantitatively check whether 
the SFH from the synthetic CMDs can explain the low measured element abundances, or if galactic 
winds with loss of metals are needed. 

Subject headings: galaxies: dwarf — galaxies: individual (I Zw 18) — galaxies: irregular — galaxies: 
resolved stellar populations — galaxies: starburst 



1. INTRODUCTION 

The blue compact dwarf (BCD) galaxy I Zw 18 is one 
of the most intriguing objects in the local Universe and 
has fas cinated genera tions of astronomers since its dis- 
covery (IZwickvi lP966l). With a metallici t y between 1/30 
and 1 /50 Z m (iSearle fc Sargent! Il972t iLeaueux et al.l 
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Skillman fc Kennicuttl P99S 
IStasihska fe Leithererl 1996; 
llzotov fc Thuanl 119981 ) it holds 
the record of the second lowest metallicity and lowest 
helium content m easured in a star-forming galaxy 
(jlzotov et al.l I2009D . The dynamical mass of I Zw 18 
measured at a radius o f 10"-12" is ~ 2-3 x 10 s M w 
(iPetrosian et all Il997t Ivan Zee et all 119981: iLelli et al.l 
|2012| ). A large amount of gas, corresponding to ~ 70% 
of the total mass, is detected all around the system, but 
only 10 7 M Q of HI is associated with the optical part of 
the g alaxy (jLequeux fc Viallefondl 119801 : Ivan Zee et al.l 
1998). The HI associated to the starburst region forms 
a compact rapidly rotating disk, and the steep rise of 
the rotation curve in the inner parts indicates that there 
is a strong central concentration of mass (jLelli et al.l 
2012). HI observations have revealed a neutral hydrogen 
bridge connecting I Zw 18 with a faint companion galaxy 
(the so called C component or secondary body), which 
has also been reso l ved into stars with HST imaging 
(IDufour et al.lll996l ; llzotov fe Thuanl [200l lAloisi et alj 
120071 ). and a ~13.5 kpc HI tail extending to the south 
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of I Zw 18 main body (jLelli et al.ll2012fl . 

I Zw 18 sho ws very blue colors, U — B = —0.88 and 
B-V = -0.03 (jvan Zee et al.lll998l) . suggesting the pres- 
ence of a very young stellar population, with a current 
star formati on rate (SFR) much h igher than the past 
mean value (jSearle fc Sargent] 1 1972f ). All these observa- 
tional characteristics make I Zw 18 resemble a primeval 
galaxy in the nearby Universe. In fact, soon after its dis- 
covery, the question arose whether I Zw 18 is so metal 
poor because a) it started forming stars only recently, so 
that they haven't had much time to pollute metals in the 
interstellar medium, or b) because its star formation (SF) 
activity, although occurring over a long period of time, 
has proceeded at a rate too low for an efficient chemical 
enrichment, or c) because strong galactic winds have re- 
moved from the system most of the metals. The nature 
of I Zw 18 has important cosmological implications. If in- 
deed some BCDs turned out to be young galaxies, their 
existence would support the view that SF in low-mass 
systems has been i nhibited till the present epoch (e.g., 
iBabul fc Fl ees 1992). On the other hand, the lack of such 
primordial systems would provide strong constraints on 
chemical evolution and hydrodynamical models of metal- 
poor galaxies, e.g., on their SF regime (continuous or 
bursting?) or the onset of a galactic wind. 

With the advent of the Hubble Space Telescope (HST) , 
it has been possible to resolve the individual stars in 
I Zw 18 and thus to charact erize its evolutionary s tatus . 
Fro m HST/WFPC2 dat a, iHunter fc Thronsonl (fl995h 
and IDufour et al.1 (|1996h argued for a continuous SF 
over the last 30-50 Myr. Thanks t o a more sophisti - 
cated treatment of the same dataset, lAloisi et al.l (jl999f ) 
were able to go deeper and to detect asymptotic gi- 
ant branch (AGB) stars with ages of at least several 
hundreds Myr. These results were subsequently con- 
firmed by lOstlinl (pOOOl ) through deep HST/NICMOS 
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Figure 1. F606W (broad V) ACS/WFC image showing I Zw 18' 
s main (south-east) and secondary (north west) bodies. 



imaging. Later, llzotov fc Thuanl (|2004D failed to detect 
red giant branch (RGB) stars with the Advanced Cam- 
era for Surveys (ACS) on board of HST and concluded 
that the galaxy is at most 500 Myr o l d. Ho wever, both 
iMomanv et"aLl (|2005l) and iTosi et ail ([20071) . from inde- 
pendent reanalyses of the same ACS data set, suggested 
that I Zw 18 should be older than at least 1-2 Gyr, since 
it does appear to contain also RGB stars. 

In order to shed light on the situation, we acquired 
in 2006 new time-series HST/ ACS photometry to study 
Cepheid stars in I Zw 18 and pin down its distance 
(GO 10586, PI Aloisi). By combining the new data with 
archival ones, we both identified the RGB ti p (TRGB) at 
Jn = 27.27 ±0.14 mag (D = 18.2±1.5 Mpc, lAloisi et al.l 
2007), and detected for the first time a few Cepheids 
whose light curves allowed us to independently derive the 
distance to ~ 1 Mpc of accuracy (D= 19 ± 0.9 Mpc, 
iFiorentino et al.l 120101: iMarconi et al.ll2010D . The detec- 
tion of RGB stars in I Zw 18 implies that it has started 
forming stars at least ~1 Gyr ago, and possibly at epochs 
as old as a Hubble time, ruling out the possibility that it 
is a truly primordial galaxy formed recently in the local 
Universe. In this paper we use synthetic CMDs to re- 
construct the entire star formation history (SFH) within 
the reachable look-back time and to quantify the stellar 
mass formed at old, intermediate-age, and young epochs. 
A qualitative study of the spatial distribution and evo- 
lutionary properties of I Zw 18' s resolved sta r s has been 
already presented bv lContreras Ramos et al.l ([201 ll ). 

2. OBSERVATIONS AND DATA REDUCTION 

Observations and data reduction have been exten- 
sively described in pre v ious papers (lAloisi et al. [20071 : 
IFiorentino et al.l 120101: iContreras Ramos et al.l I2011D . 
Here we only recall the main aspects. The data were 
collected between October 2005 and January 2006 with 
the ACS/WFC in 13 different epochs in F606W (-broad 
V), and 12 different epochs in F814W (~I), for total in- 
tegration times of « 27,700 s and m 26,200 s. The single 
exposures per epoch per filter were reprocessed with the 
most up-to-date version of the ACS calibration pipeline 
(CALACS). Then, for each filter, we co-added all the 
exposures into a singl e image using the softw are pack- 
age MULTIDRIZZLE (|Koekemoer et al.ll2002h . in order 
to obtain 2 deep images in F606W and F814W, respec- 
tively. A portion of the F606W deep image showing 
I Zw 18' s main and secondary bodies is shown in Fig. [TJ 
The pixel size of the final resampled drizzled images is 



0.035 "(0.7 times the original ACS/WFC pixel size). We 
also retrieved and combined archival ACS/WFC data 
in F555W (~V) and F814W (GO program 9400, PI 
Thuan), taken in 2003 over a period of 11 days, into 
two deep master images with total integration times of 
« 43,500 s and « 24,300 s, respectively. 

Photometry was performed with DAOPHOT (|Stetsonl 
1987) in the IRAF environment for both the archival 
and proprietary datasets. The instrumental magnitudes 
were estimated via a PSF-fitting technique. The PSF was 
created choosing the most isolated and clean stars in the 
vicinity of I Zw 18's main and secondary bodies, and was 
modeled with an analytic Moffat function plus additive 
corrections derived from the residuals of the fit to the 
PSF stars. In order to push the photometry as deep as 
possible, stars were detected in a sum of the V and I im- 
ages in each dataset, choosing a detection treshold of 2.4 
times the local background level. This procedure allowed 
us to recover faint objects with very blue or red colors. 
Aperture photometry with PHOT, and then PSF-fitting 
photometry w ith the DAOPHOT/ ALLSTAR package 
(|Stetsonl lT987). were performed separately on the V and 
I images at the position of the objects detected in the 
sum image. The finding procedure, the aperture pho- 
tometry, and the PSF-fitting photometry were then re- 
iterated on the subtracted V and I images. This al- 
lowed us to recover additional faint stars showing up 
only after the subtraction of their brighter neighbours. 
For each dataset, the HST magnitudes were calibrated 
into the Johnson-Cous i ns sys tem using the transforma- 
tions in iSirianni et al.l (|2005| ). The CMD presented in 
lAloisi et al.1 (|2007l ) was obtained by cross-correlating the 
photometry in the two datasets requiring good matches 
in position (within one pixel) and magnitude (within 
3.5cr), and a Daophot sharpness parameter < 0.5. The 
final V and I magnitudes were obtained averaging the 
magnitudes in the two datasets. In this paper, we use 
this matched catalog with «2000 sources to derive the 
SFH of I Zw 18. 

3. COLOR-MAGNITUDE DIAGRAMS 

To account for the highly variable crowding within 
I Zw 18, we selected five regions (A, B, C in the main 
body, D and E in the secondary body) . The correspond- 
ing I, V— I color-magnitude diagrams (CMDs) are shown 
in Fig.[2j In the main body, region A is the most crowded 
one, and is embedded into a shell of ionized gas, while re- 
gion C, at the outskirts, is the least crowded one. In the 
secondary body, region E is more crowded than region 
D. We detect 111, 481, 590, and 548, and 265 sources 
in regions A, B, C, D, and E, over an area of ~0.06, 
0.66, 4.83, 0.17, and 0.70 kpc 2 , respectively. However, 
in region C the majority of stars are concentrated into 
a smaller area of 0.76 kpc 2 , and thus we will refer to 
this value hereafter to compute the specific star forma- 
tion rate ( see Table ffl). The 8. 7 -day period classical 
Cepheid of IFiorentino et all (|2010[) is located in region 
C, the o ther two Ultra long-pe riod Cepheids (P—130 d, 
see also IFiorentino et al.l I2012D are found in regions C 
and A, while the two long period variables (P>100 d) 

5 IRAF is distributed by the National Optical Astronomy Ob- 
servatories, which are operated by AURA, Inc., under cooperative 
agreement with the National Science Foundation 
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are in regions B and C. 

On the CMDs in Fig. [2l we have overplotted the 
Padova 94 t F agotto et all 19941 ) stellar tracks for different 
masses to show the mass ranges (and thus the look-back 
times) sampled in the different regions. In region A, a 
concentration of bright stars between the 20 and 30 M@ 
tracks, with ages 7-10 Myr, and a few stars with mass 
as low as ~ 12 M Q (^20 Myr old) are detected. Lower 
mass stars are lost because of the severe photometric in- 
completeness due to the extreme crowding. In region B, 
the look-back time sampled is longer than in region A, 
with stellar masses as low as 2 M Q (~ 1 Gyr old) present 
on the CMD. However, only in the least crowded region 
C red giant branch (RGB) stars with masses < 2M Q and 
ages > 1 Gyr are detected. This also holds for region D 
in the secondary body, while in the more crowded region 
E only a few potential RGB stars are observed. From 
the comparison of the CMDs with the theoretical tracks, 
we also notice the clear presence of an age gradient mov- 
ing from the most crowded regions to the least crowded 
ones. In fact, the bright part of the blue plume (I<26 
mag) is significantly more populated in regions A, B and 
E than in regions C and D, and reaches brighter mag- 
nitudes. This implies more massive stars and younger 
star formation. For a quantitative derivation of the star 
formation history (SFH) in the different regions we refer 
to Section [5l 

4. ARTIFICIAL STAR EXPERIMENTS 

To evaluate the role of incompleteness and blending 
in the data, we performed artificial star experiments on 
the origin al frames, following t he same procedure de- 
scribed in lAnnibali et al.l (|2008). These tests serve to 
probe observational effects associated with the data re- 
duction process, such as the accuracy of the photometric 
measurements, the crowding conditions, and the ability 
of the PSF-fitting procedure in resolving partially over- 
lapped sources. We performed the tests for the individual 
F555W, F606W, and for the two F814W frames, accord- 
ing to the following procedure. We divided the frames 
into grids of cells of chosen width (50 pixels) and ran- 
domly added one artificial star per cell at each run. This 
procedure prevents the artificial stars to interfere with 
each other, and avoids to bias the experiments toward 
an artificial crowding not present in the original frames. 
The position of the grid is randomly changed at each 
run, and after a large number of experiments the stars 
are uniformly distributed over the frame. In each filter, 
we assign to the artificial star a random input magnitude 
between ml and m2, with ml « 3 mag brighter than the 
brightest star in the CMD, and m2 « 3 mag fainter than 
the faintest star in the CMD. At each run, the frame is re- 
reduced following the same procedure as for the real data. 
To account for the fact that in the real data the source 
detection was performed on a sum of the V+I frames, we 
adopted a smaller treshold than the value of 2.4 adopted 
for the real data. The new treshold is derived performing 
tests on the original frames, and is such that the number 
of detections in a single frame equals the number of de- 
tections in the sum of the V+I frames with a treshold of 
2.4. The output photometric catalog is cross-correlated 
with a sum of the original photometric catalog of real 
stars and the list of the artificial stars added into the 
frame. This prevents cross-correlation of artificial stars 



in the input list with real stars recovered in the output 
photometric catalog. We simulated ^500,000 stars for 
each filter. Stars with input-output magnitude > 0.75 
were considered lost, because such a difference implies 
that they fell on a real star of their same luminosity or 
brighter. For each magnitude bin, the completeness of 
our photometry was computed as the ratio of the num- 
ber of recovered artificial stars over the number of added 
ones. 

As an example, we show in Fig. [3] the distribution of 
the artificial star output minus input magnitude (Am) in 
F814W for the Aloisi dataset. We show the different dis- 
tributions relative to regions A, B, C, D and E selected in 
I Zw 18, in order to characterize the photometric error as 
a function of magnitude and crowding. The systematic 
deviation from of the mean Am is due to the increasing 
effect of blending at fainter magnitudes: faint stars tend 
to be systematically recovered with brighter magnitudes 
because they happen to overlap with objects of compara- 
ble or brighter magnitude. The photometric errors sys- 
tematically decrease from the most crowded region A to 
the least crowded regions C and D. At the same time, 
the completeness increases, as shown in Fig. 01 

5. STAR FORMATION HISTORY 

5.1. Methodology 

We derived the star formation history (SFH) of I Zw 18 
through comparison of the observed CMDs with syn- 
thetic ones. Our procedure consists of two separate steps 
as described in the following subsections: 1) the defini- 
tion of the so called basis functions from the combination 
of the adopted stellar evolution models with the pho- 
tometric properties of the examined region, and 2) the 
statistical analysis for the identification of the best SFH 
and relative uncertainties. In order to better characterize 
the uncertainties intrinsic in the statistical analysis, we 
applied two completely independent procedures: the one 
described bv lGrocholski et all (|2012[ ) (he reafter, the Bal- 
timore procedure) and that described bv lCignoni fc Tosll 
(2010) (hereafter, the Bologna procedure). The compar- 
ison between the two methods provides important in- 
formation on which features of the SFH are robust and 
which are instead artifacts due to the adopted minimiza- 
tion approach. On the other hand, quantifying the sys- 
tematics due to the choice of a particular set of stellar 
models goes beyond the scope of this paper. 

5.1.1. Basis Functions 

First, we created a grid of synthetic CMDs for episodes 
with log(age) from 6 to 10, and duration Slog(age) = 
0.25 (hereafter, basis functions). The choice of logarith- 
mic steps allows us to account for the increasingly lower 
time resolution in the SFH at larger look-back time. The 
basis functions are created fo l lowin g the approach ini- 
tially described in iTosi et all (119911), and a d optin g the 
last version of the code by lAngeretti eFall (|2005l ). In 
brief, the synthetic CMDs are produced via Monte Carlo 
extractions of (mass, age) pairs for an assumed Initial 
Mass Function (IMF), adopting a constant SF from the 
starting to the ending epoch of each basis func tion. Here 
we adopt a Salpeter' s IMF ()Salpeterl [1955) and zero 
binary fraction. Each star is placed in the theoretical 
(logL/L , logTe//) plane via interpolation on a chosen 
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Figure 2. Upper panels: spatial distribution of the stars in regions A, B, C, D and E selected in I Zw 18. Bottom panels: corresponding 
CMDs. Superimposed are the Padova 94 stellar tracks for the masses: 60, 30, 20, 12, 7, 5, 2, 1, 0.8, and 0.6 Mq. For regions A and B, the 
smallest plotted masses are 20 and 2 Mq, respectively. 



set of stellar evolutio nary tracks (here, the Padova tracks 
with Z = 0.0004, iFa'gotto et all 119941 ) . Luminosity 
and effective temperature are transformed into the ACS 
Vegamag F555W, F606W and F814W filters using the 
lOriglia fc Leithererl (|2000f ) code. Absolute magnitudes 
are then converted into apparent ones applying a fore- 
ground reddening of E(B-V) = 0.032 from the NED and 
a distance modulu s of (m — Mo) = 31.3 mag derived by 
lAloisi et~al"1 ((2007). At this point we apply a complete- 
ness test in order to determine whether to retain or to 
reject the synthetic star, based on the results of the arti- 
ficial star tests performed in all four images (F606W and 
F814W for the Aloisi dataset, and F555W and F814W 
for the Thuan dataset, see Section @|. Since we are us- 
ing the CMD obtained from the cross-correlation of four 
photometric catalogs (see Section \S§ , we require in the 
simulations that the synthetic stars pass the test in all 
four photometric bands. Photometric errors are assigned 
on the basis of the distribution of the output-minus-input 
magnitudes of the artificial stars (Fig. [3]). These errors 
take into account the various instrumental and observa- 
tional effects, as well as systematic uncertainties due to 
crowding (i.e., blend of fainter objects into an apparent 
brighter one). Four magnitudes (F555W, F606W, plus 2 
in F814W) are associated to each synthetic star. The ex- 
traction of (mass, age) pairs is stopped when the number 
of stars populating the synthetic CMD equals a certain 



fixed number. This number was chosen to be the same 
for all the basis functions and to guarantee that all the 
phases potentially sampled in the synthetic CMD (given 
the completeness function) turn out to be well populated. 
For a direct comparison with the observed CMDs in Sec- 
tion O the ACS Vegamag magnitudes are transformed 
into the Johnson- Cousins system and are averaged fol- 
lowing the same procedure as for the real data. The same 
basis functions are implemented into the Baltimore and 
the Bologna procedures. 

5.1.2. Statistical analysis with the Baltimore approach 

The SFH was derived through a statistical approach 
using the code SFHMATRIX develop ed by R. van der 
Mare l in Baltimore and described in ()Grocholski et all 
I20l2j) . In brief, the SFH is inferred finding the weighted 
combination of basis functions that best reproduce the 
observed CMD in a x 2 sense. To maximize the likeli- 
hood, the code considers the density of points on the 
observed/synthetic CMD, i.e. the Hess diagram, rather 
than the CMD itself. The best fit is found by solving 
a non-negative least-squares matrix problem. In SFH- 
MATRIX, the errors on the SFH are calculated creating 
many realizations of pseudo data sets with properties 
similar to the real data, and analyzing these pseudo data 
in the same way as for the real data. Then, the rms scat- 
ter in the SFH for a given age is the error bar. We create 
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Figure 3. Output minus input magnitude versus input magnitude in F814W for the artificial stars simulated in the Aloisi dataset. The 
distribution is shown for the different regions (A, B, C, D, E) selected in the 1 Zw 18' s field. The solid curves indicate the mean of the 
Amag distribution (central line) and the ±lcr standard deviations. For regions B, C, D and E, we give the a value at F814W=28. 



the pseudo-data by drawing many Monte-Carlo realiza- 
tions from the best-fit SFH that was inferred from the 
data. 

5.1.3. Statistical analysis with the Bologna approach 

We also derived the SFH using the Bologna code (fo r 
details and application see Cien oni et all l201lL \20m . 
which has many features in common with the Baltimore 
code, but also incorporates some differences in the search 
of the best synthetic CMD. Similarly to the Baltimore 
algorithm, the SFH is parametrized as a linear combi- 
nation of quasi simple stellar populations with variable 
duration (the basis functions), generated from a specific 
set of stellar models with given metallicity. The code 
searches for the linear combination of the basis functions 
that best matches the observed CMD. The comparison 
between the observed and the synthetic CMDs is per- 
formed through a \ 2 minimization in strategic regions 
of the CMD. As the Baltimore code, the Bologna code 
considers the Hess diagram rather than the CMD itself 



to perform the minimization. The sizes and distribu- 
tion of the bins within the color-magnitude plane can be 
variable, and are chosen considering some aspects such 
as the particular CMD shape, the number statistics, the 
presence of phases where the stellar evolution models are 
more uncertain, and so on. Differently from the Balti- 
more code, the Bologna code employs for the minimiza- 
tion a downhill simplex "amoeba" algorithm which is re- 
started from thousands of initial random positions. A 
bootstrap method is used to assess the statistical uncer- 
tainty around the best-fit parameters. The minimization 
process is repeated for each bootstrapped data set, and 
the final error bars represent the mean deviation using 
1000 bootstraps. 

5.2. CMD fit 

5.2.1. Region A 

Fig. shows the comparison between the observed 
CMD selected in the most crowded region A (left) and 
two synthetic CMD realizations drawn from the best- 



6 



E 







1 

0.8 


- 




\ ---reg C - 

\. 1 reg D - 

v \ \ - - reg E - 

\ \ \ - 


0.6 


: 




Ml - 
\\ 


0.4 






\ \ 


0.2 






U 

w 









V v 




i 





20 25 

F814W (Aloisi) 



Figure 4. Completeness factor (1=100% complete) as a function 
of the F814W magnitude for the Aloisi dataset. The different lines 
correspond to the different regions A, B, C, D, and E selected in 
I Zw 18, as indicated by the labels. 

fit SFH with the Baltimore code (middle) and with the 
Bologna code (right). In both cases, when searching for 
the best-fit solution, we excluded the region below the 
60% completeness limit from the \ 2 minimization. The 
reason is that crowding is extreme in region A, implying 
that the derived completeness curve may become highly 
uncertain below this limit. When running the Baltimore 
code, we adopted a constant pixel size of 0.5 in both color 
and magnitude. Since the stellar tracks fail to reach the 
reddest supergiants in the observed CMD, we decided to 
treat the box at 0.8 < V-I < 2, 21.5 < I < 23 as a single 
pixel. On the other hand, in the Bologna procedure, the 
CMD was divided into a regular grid of color-magnitude 
cells with A(V - I) = 0.25 mag and AI = 1 mag. The 
choice of the pixel size, which is driven by the need of 
minimizing statistical fluctuations without loosing the in- 
formation in the CMD, remains somewhat subjective. 
Thus, it is interesting to see how the results change ac- 
cording to the different choices of the grid. Both the 
Baltimore and the Bologna codes provided reduced % 2 
values less than 1 for the best-fit solutions. The absolute 
X 2 of the fit gives a measure of how well the data are re- 
produced, but it is often difficult to interpret in practice. 
For example, the effective number of degrees of freedom 
is difficult to estimate when many Hess diagram pixels 
are empty (such as in region A), and do not effectively 
contribute to constraining the SFH. We therefore prefer 
to assess the quality of fit in this paper through other 
metrics, including data-model comparisons of luminosity 
functions. 

Fig.[5]shows that both the Baltimore and the Bologna 
solutions well reproduce the overall morphology of the 
observed CMD; however, an inspection of the luminos- 
ity functions (LFs) in Fig. |6] reveals some discrepancies 
between observations and simulations. More specifically, 
the Baltimore solution tends to underpredict the blue 
(V - I < 0.6) counts at 22 < I < 22.5, and the red 
(V — I > 0.6) counts at 22 < I < 23; these stars are 
blue and red supergiants with masses higher than ~ 15 
Mq and ages less than ~10 Myr in their post-main se- 
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Figure 5. From left to right: observed CMD for region A, best- 
fit synthetic CMD with the Baltimore code, best-fit CMD with 
the Bologna code. The best-fit CMDs were obtained considering 
only the region above the 60% completeness line, indicated by the 
dotted curve. The box includes red supergiants with age 7-10 Myr. 
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Figure 6. Luminosity functions in the blue (left) and in the red 
(right) for the observed CMD in region A (points) and the best-fit 
CMDs obtained with the Baltimore (solid line) and the Bologna 
(dotted line) codes. 

quence (post-MS) phase (we will ignore the shortage of 
simulated stars at I > 25.5 mag, since this region was 
excluded from the fit and the completeness correction 
may be uncertain). Also the Bologna solution tends to 
underpredict the number of red supergiants, but pro- 
vides a better match to the bright blue supergiants at 
22 < I < 22.5. However, this occurs at the expenses 
of an excess of fainter blue counts at 24 < / < 25. 
Concluding, we notice that both solutions tend to un- 
der reproduce the correct ratio of brighter over fainter 
massive post-MS stars. As shown in Fig. [71 where we 
have coded the mass in the simulated CMD, this prob- 
lem can not be solved with a flatter IMF, because the 
same masses populate at the same time both bright and 
faint regions in the CMD. Possible explanations are: a) 
the completeness behaviour as a function of magnitude is 
not properly taken into account; indeed the severe crowd- 
ing in region A may prevent a robust determination of 
the completeness function through artificial star tests; 
b) some of the brightest stars in the CMD are in fact 
blends of two or more stars or unresolved clusters rather 
than individual stars; however, as already discussed in 



7 



21 



22 - 



1 1 1 1 1 1 1 1 1 1 


1 1 1 

o<* 


1 1 1 1 1 1 1 1 1 1 
MODEL - BALTIMORE: 
REG A 


_ x 

* x 

o 


X x 

X 




Oft, x 






- J3T it 


• 

• 


•10-20 M '. 

x 20-30 M - 


^x* 

X xx x 

X X 

X 




030-40 M '. 
*40-50 M - 


"I 1 1 1 1 1 1 ? 1 Ift 


1 1 1 


A >50 M '. 
■ i i 1 i i i i 1 i~ 


-1 


1 

V-I 


2 3 



Figure 7. Best-fit CMD with the Baltimore code for region A. 
Different symbols correspond to different masses, as indicated by 
the labels. 
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Figure 8. Best-fit SFH for region A with the Baltimore code 
(left panels) and with the Bologna code (right panels). From top 
to bottom: stellar mass formed and SFR in the different age bins. 



IContreras Ramos et al.1 (|2011| ). the colors of the reddest 
(V — I ^ 1) bright objects in I Zw 18 are incompatible 
with the cluster hypothesis; we can also exclude that the 
large number of bright red stars observed in the CMD 
is due to reddening within I Zw 18; in fact, the largest 
reddening value estimated bv lCannon et al.l ([2002D in the 
NW component from the Ha/H/3 flux ration amounts to 
just E(V — J) ~ 0.1; c) the evolutionary timescales in 
the brightest post-MS phase for massive stars are un- 
derestimated. The same discrepancy persists if we run 
some simul ations with the new Padova PARSEC stel- 
lar models (|Bressan et al.ll201~2l and Bressan et al., pri- 
vate communication) created with a thorough revision of 
the major input physics, a new treatment of the opaci- 
ties consistent with the chemical composition, and a new 
treatment of the mass loss. It is possible that models in- 
cluding additional input physics, such as stellar rotation 
(|Mevnet &; Ma edcr 2002), could provide a better match 
to the brightest objects in the CMD, but these models 
are not available at the moment for a large set of masses 
at I Zw 18's metallicity. 



The best-fit SFHs obtained with the two different pro- 
cedures are shown in Fig. [8] For both the Baltimore and 
the Bologna code, we show the stellar mass formed and 
the SFR in the different age intervals. The results are 
also summarized in Tables [1] and [21 The values are nor- 
malized to a Salpeter IMF in the mass interval 0.1 to 
120 Mq. The symbol > implies that the value is a lower 
limit. This is because the look-back time may not be 
large enough to probe the entire age-range in some re- 
gions. The look-back time in region A is ~20 Myr. Both 
solutions show a strong SF event occurred 6-10 Myr ago 
with a rate as high w 1.5M p/yr, comparable to the high 
rate derived in NGC 1569 (|Greggio et al.lll998[ ). A mass 
of r; 6 x 1O 6 M , similar to that of the largest globular 
clusters ever discovered (e.g. uj Cen in the Milky Way 
and Gl in M 31) was formed during this episode. If the 
SFR is divided by the area sampled by Region A, we ob- 
tain ~ 2.5 x W MQyr pc , e.g. just a few factors 
lower than the rate derived in the strong star forming 
region 30 Doradus in the LMC. From our best- fit SFH 
we can not exclude the presence of an even more recent 
burst, or ongoing SF. 

5.2.2. Region B 

The comparison between the observed CMD in region 
B and two synthetic CMD realizations from the best-fit 
SFH with the Baltimore and the Bologna codes is shown 
in Fig. [§J We excluded from the fit the area below the 
dotted line where the completeness is below 40%. Given 
the higher density of stars on the CMD than in region A, 
we adopted here a smaller pixel size of 0.25 in both color 
and magnitude with the Baltimore code, while we chose 
a pixel size of 0.25 in color and 0.5 in magnitude to run 
the Bologna code. When running the Baltimore code, 
we treated the box at 0.8 < V - I < 2, 21.5 < I < 23 
as a single pixel (as for region A) , and excluded the box 
at 0.8 < V - I < 2, 24 < I < 27, where TP-AGB 
stars are located, from the fit. The reason for this is 
that the TP-AGB phase is not included in the Padova 94 
tracks that we use i n this paper. Altho ugh more recent 
Padova models (see IGirardi et al.|[20Toh provide a much 
improved treatment of the TP-AGB stars, this is still a 
complex and highly uncertain phase of stellar evolution. 
Moreover, a significant but unknown fr action of TP-AGB 
stars may be obscured by dust shells ()Bover et al.ll2009l) 
making it difficult to use such stars to constrain the SFH. 
On the other hand, the Bologna code was run considering 
the entire CMD down to the 40% completeness limit. We 
purposely adopted a different approach for the TP-AGB 
region in the two codes to be able to assess the impact 
of this region on the final inferred SFH. This is valid for 
all the other regions as well. 

As shown in Fig.[5J both solutions reproduce quite well 
the overall morphology of the observed CMD. The agree- 
ment is reasonable also when comparing the observed 
and simulated LFs in Fig. [TO] The main discrepancies 
that we notice are the following: 1) both solutions tend 
to under-reproduce the number of bright blue stars at 
/ < 24 mag and bright red stars 22 < / < 23 (we noticed 
an analogous problem in region A); 2) the Baltimore so- 
lution tends to under reproduce the faint red counts at 
I > 27.5 mag, while the Bologna solution over predicts 
the blue counts at / > 27.5 mag; however, we recall that 
these regions fall below the 40% completeness limit and 
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Figure 9. From left to right: observed CMD for region B, best- 
fit synthetic CMD with the Baltimore code, best-fit CMD with 
the Bologna code. The best-fit CMDs were obtained considering 
only the region above the 40% completeness line, indicated by the 
dotted curve. The upper box includes red supergiants with age 
7-10 Myr. The lower box is where TP-AGB stars are located. 
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Figure 10. Luminosity functions in the blue (left) and in the red 
(right) for the observed CMD in region B (points) and the best-fit 
CMDs obtained with the Baltimore (solid line) and the Bologna 
(dotted line) odes. 





: reg B - BALTIMORE 


1 1. 


T 10" 






3 






g 5xl0 5 








3 















, h 



log(age[yr]) 



log(age[yr]) 




log(age[yr]) 



log(age[yr]) 



Figure 11. Best-fit SFH for region B with the Baltimore code 
(left panels) and with the Bologna code (right panels). From top 
to bottom: stellar mass formed and SFR in the different age bins. 



were excluded from the fit when searching for the best-fit 
solution. 

The best-fit SFHs obtained with the two different pro- 
cedures are shown in Fig.QTJ In both solutions, the SFR 
decreases with increasing look-back time (—300 Myr is 
the maximum look-back time reached in this region) ; the 
highest peak of SF occurred at a rate of ~ 0.04M Q /yr 
during the last 10-20 Myr (depending on the solution, 
see Table [IJ . The total mass formed in region B over the 
last -300 Myr is w 2 - 3 x 10 6 M Q , depending on the 
solution (see Table [2]) . 



5.2.3. Region C 

As shown in Fig. [2] the CMD of region C is deep 
enough to reach RGB stars and thus to sample the SF at 
epochs older than —1 Gyr. The comparison between the 
observed CMD and two synthetic CMD realizations from 
the best-fit SFH with the Baltimore and the Bologna 
codes is shown in Fig. [T2J Here we were less conservative 
than in regions A and B, and performed the % 2 minimiza- 
tion using the CMD down to the 20% completeness limit, 
in order to reach faint stars in the RGB phase (I > 27.5 
mag). When running the Baltimore code, we adopted a 
pixel size of 0.25 in magnitude and color and, as for region 
B, we excluded from the fit the box at 0.8 < V — I < 2. 
24 < I < 27 where TP-AGB stars are located. On the 
other hand, the Bologna code was run adopting a pixel 
size of 0.5x0.5 and 0.25x0.25 in the magnitude ranges 
22 < I < 25 and 25 < I < 29, respectively, and without 
the exclusion of CMD regions down to the 20% complete- 
ness limit. Both solutions provide a good agreement with 
the overall CMD morphology in Fig.[T2land with the LFs 
in Fig. [13] The Bologna solution tends to exceed the 
red counts at 28 < I < 28.5, while the Baltimore solu- 
tion tends to under predict the red counts at the faintest 
(/ > 28.5) magnitudes; we recall that this region, where 
the completeness is very uncertain, was excluded from 
the fit. In general, both solutions provide a good agree- 
ment with the observed counts at the RGB tip (I — 27.5 
mag); here the completeness is ~70%, high enough to 
allow a reliable derivation of the old (> 1 Gyr) SFH. 

The best-fit SFHs obtained with the Baltimore and 
the Bologna procedures are shown in Fig. [T5J Accord- 
ing to both solutions, the average SFR over the last 100 
Myr was «10 times as large as the average rate over the 
1 — 10 Gyr interval (see Table Q] for details). However, 
the lack of resolution at epochs older than 1 Gyr (at a 
metallicity of Z = 0.0004 the RGB color is highly de- 
generate with age and our photometric errors are large) 
does not allow us to establish if short bursts at a much 
higher rate occurred in I Zw 18. Indeed, despite the 
fact that we provide the SFR in four different age inter- 
vals from log(a<7e) =9 to 10, the only meaningful quan- 
tity at these ages is the total stellar mass formed, i.e 
?» (0.9 — 1.4) x 1O 6 M , depending on the solution (see 
Table [5]) . This is also evident from the comparison of 
the SFHs with the Bologna and the Baltimore solutions, 
which differ somewhat from each other prior to 1 Gyr 
ago, but still within the uncertainties. 

Given that the RGB tip is so close to the completeness 
limit, we also tested the impact of the distance modulus 
on the final SFH. We fo und that, w it hin t he distance 
uncertainties derived by lAloisi et all (|2007t ) (18.2±1.5 
Mpc), there are no significant effects on the final SFH. 
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Figure 12. From left to right: observed CMD for region C, best- 
fit synthetic CMD with the Baltimore code, best-fit CMD with 
the Bologna code. The best-fit CMDs were obtained considering 
only the region above the 20% completeness line, indicated by the 
dotted curve. The box is where TP-AGB stars arc located. 
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Figure 15. From left to right: observed CMD for region D, best- 
fit synthetic CMD with the Baltimore code, best-fit CMD with 
the Bologna code. The best-fit CMDs were obtained considering 
only the region above the 20% completeness line, indicated by the 
dotted curve. The box is where TP-AGB stars arc located. 
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Figure 13. Luminosity functions in the blue (left) and in the red 
(right) for the observed CMD in region C (points) and the best-fit 
CMDs obtained with the Baltimore (solid line) and the Bologna 
(dotted line) odes. 
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Figure 14. Best-fit SFH for region C with the Baltimore code 
(left panels) and with the Bologna code (right panels). From top 
to bottom: stellar mass formed and SFR in the different age bins. 



5.2.4. Region D 

The CMD of region D in the secondary body ap- 
pears quite similar to that of region C in the main body 
(Fig. [2]). The photometry is deep enough to potentially 
reach RGB stars sampling the SF at epochs older than 
~1 Gyr. The comparison between the observed CMD 
and two synthetic CMD realizations from the best-fit 
SFH with the Baltimore and the Bologna codes is shown 
in Fig. [15] As for region C, the x 2 minimization was done 
using the CMD down to the 20% completeness limit, in 
order to reach faint stars in the RGB phase. When run- 
ning the Baltimore code, a pixel size of 0.25 in both color 
and magnitude was adopted and the region of TP-AGB 
stars at 0.8 < V - I < 2, 24 < I < 27 was excluded 
from the fit. On the other hand, with the Bologna code 
we adopted a pixel size of 0.25x0.25 and 0.5x0.5 in the 
magnitude ranges 22 < / < 24 and 24 < I < 29, re- 
spectively, and used all the CMD down to the 20% com- 
pleteness limit for the % 2 minimization. Both solutions 
reproduce quite well the overall CMD morphology. The 
agreement with the LFs in Fig. \W\ is quite good too, 
with the only exception of the region at 27.5 < I < 28.5, 
V — I > 0.6, where the Baltimore solution tends to over 
predict the counts. 

The best-fit SFHs obtained with the two different pro- 
cedures are in very good agreement with each other 
(Fig. [17]) ■ In both solutions, the highest peak of SF oc- 
curred ~ 250 Myr ago at a rate of ~ 5 x 10~ 3 Mo /yr, 
and the SFR has been gradually decreasing since then. 
The Baltimore solution is compatible with SF during the 
1 — 10 Gyr age interval, with an upper limit of 10 6 A/q 
for the stellar mass (see Figure IT7|). while the Bologna 
solution seems to exclude the presence of significant SF 
at these epochs. This implies that the fraction of the 
mass in old stars in region D is poorly constrained, and 
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Figure 16. Luminosity functions in the blue (left) and in the red 
(right) for the observed CMD in region D (points) and the best-fit 
CMDs obtained with the Baltimore (solid line) and the Bologna 
(dotted line) odes. 
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Figure 17. Best-fit SFH for region D with the Baltimore code 
(left panels) and with the Bologna code (right panels). From top 
to bottom: stellar mass formed and SFR in the different age bins. 

might be anywhere between and 40%. 

5.2.5. Region E 

The comparison between the observed CMD in the 
most crowded region E of the secondary body and two 
synthetic CMD realizations from the best-fit SFH with 
the Baltimore and the Bologna codes is shown in Fig. [TBI 
Due to the higher crowding than in region D, it is more 
difficult here to sample RGB stars. In fact, the com- 
pleteness is as low as ^30% at the RGB tip. The fit 
was done considering the CMD down to the 20% com- 
pleteness limit. When running the Baltimore code, we 
excluded the region of the TP-AGB stars, and adopted 
a pixel size of 0.25 in both color and magnitude. On the 
other hand, when running the Bologna code, we consid- 
ered all the CMD down to the 20% completeness limit 
and adopted a pixel size of 0.5x0.5 and 0.25x0.25 in the 
magnitude ranges I < 26 mag and I > 26 mag, respec- 
tively. 

From Figs. IT51 and HT)1 we notice that the Baltimore so- 
lution tends to under predict the blue (V — I < 0.6) 
counts at I > 28 mag while overproducing the red 
(V - I > 0.6) counts at 24.5 < I < 27.5; since this 
is the region of TP-AGB stars, that was excluded from 
the fit, the poor agreement between simulations and ob- 
servations should not surprise. On the other hand, the 
Bologna solution seems to better reproduce both the 
overall CMD morphology and the LFs. Despite these 
differences, the best-fit SFHs obtained with the two pro- 
cedures (Fig. I20p are in very good agreement with each 
other. The highest peak of SF occurred 30-10 Myr ago 
at a rate of ~ 7 x lO _3 M /yr. This is much lower than 
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Figure 18. From left to right: observed CMD for region E, best- 
fit synthetic CMD with the Baltimore code, best-fit CMD with 
the Bologna code. The best-fit CMDs were obtained considering 
only the region above the 20% completeness line, indicated by the 
dotted curve. The lower box is where TP-AGB stars arc located. 





Figure 19. Luminosity functions in the blue (left) and in the red 
(right) for the observed CMD in region E (points) and the best-fit 
CMDs obtained with the Baltimore (solid line) and the Bologna 
(dotted line) odes. 
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Figure 20. Best-fit SFH for region E with the Baltimore code 
(left panels) and with the Bologna code (right panels). From top 
to bottom: stellar mass formed and SFR in the different age bins. 
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Average SFR (in 10 



Table 1 

" 3 M /yr) in I Zw 18 at different 
epochs 



Region a 


Area 


<10 


10-100 


0.1-1 


1-10 




[kpc 2 ] 


[Myr] 


[Myr] 


[Gyr] 


[Gyr] 


A 
i\ 


U.UD 


11U ( . U 












[1026.0] 


>[14.9] 




- 


B 


0.66 


20.0 


12.7 


>0.7 








[23.0] 


[16.2] 


>[2.5] 




C 


0.76 


2.0 


1.5 


0.8 


0.2 






[2.0] 


[2.1] 


[1.4] 


[0.1] 


MB 


1.48 


1129.0 


>38.9 


>1.5 


>0.2 






[1051.0] 


>[33.2] 


>[3.9] 


>[0.1] 


D 


0.17 


0.3 


1.3 


1.1 


0.04 






[0.4] 


[1.4] 


[1.0] 


[<0.01] 


E 


0.70 


0.5 


3.8 


0.4 


0.03 






[2.0] 


[3.0] 


[0.5] 


[<0.01] 


SB 


0.87 


0.8 


5.1 


1.5 


0.07 






[2.4] 


[4.4] 


[1.5] 


[<0.02] 


Note. 


— Values 


are given 


for the 


Baltimore solution 



and, in brackets, for the Bologna solution. The symbol > 
is for lower limits. 

a Region: MB=main body; SB=secondary body. 



Table 2 

Stellar mass (in 1O 6 M ) formed in I Zw 18 at different 
epochs 



Region 3 


<10 


10-100 


0.1-1 


>1 


Total 


[Myr] 


[Myr] 


[Gyr] 


[Gyr] 




A 


11.07 


>2.22 






>13.29 




[10.26] 


>[1.34] 






>[11.60] 


B 


0.20 


1.14 


>0.61 




>1.95 




[0.23] 


[1.46] 


>[1.68] 




>[3.37] 


C 


0.02 


0.14 


0.70 


1.45 


2.31 




[0.02] 


[0.19] 


[1.25] 


[0.87] 


[2.33] 


MB 


11.29 


>3.50 


>1.31 


>1.45 


>17.55 




[10.51] 


>[2.99] 


>[2.93] 


>[0.87] 


> [17.30] 


D 


0.003 


0.12 


0.98 


0.37 


1.47 




[0.004] 


[0.13] 


[0.87] 


[0.01] 


[1.01] 


E 


0.005 


0.34 


0.40 


0.28 


1.02 




[0.020] 


[0.27] 


[0.48] 


[0.02] 


[0.79] 


SB 


0.008 


0.46 


1.38 


0.65 


2.50 




[0.024] 


[0.40] 


[1.35] 


[0.03] 


[1.80] 



Note. — Values are given for the Baltimore solution and, 
in brackets, for the Bologna solution. The symbol > is for 
lower limits. 

a Region: MB=main body; SB=secondary body. 

what found for the most active regions of the main body. 
The Baltimore solution is compatible with the presence 
of SF at ages older than 1 Gyr (see Tables [Hand [2|). 

6. DISCUSSION 

In p revious papers (lAloisi et al.ll2007t iFiorentino et al.1 
2010; iContreras Ramos et al.l 1201 ID we have tried to 
qualitatively characterize the different stellar populations 
in I Zw 18 and have demonstrated the presence of RGB 
stars, excluding the possibility that this is a truly primor- 
dial system in the present-day Universe. In this paper 
we use HST/ACS data in combination with the method 
of synthetic CMDs to quantitatively derive the SFH in 
I Zw 18. Prev i ous si milar st udies have b e en pr esented 
bv lAloisi etldl (| 19991) an d bvlJamet et all (120101 ) . How- 
ever, the SFH derived by lAloisi et al.l "([l999) was based 



on shallower WFPC2 data and on the assumption of a 
distance of 10 Mpc, much closer than what is suggested 
by the deeper ACS CMP and by t he Cepheid study. On 
the other hand. lJamet et all (|201CO used ACS data (only 
the Thuan dataset) for the derivation of the SFH in the 
secondary body of I Zw 18, but inferred a distance of 
~27 Mpc from the CMD analysis, incompatible with the 
distance of ^18 Mpc obtained from the Cepheid light 
curves. In this paper we present for the first time the 
SFH of I Zw 18 based on a robust distance assumption. 

In order to account for the different crowding condi- 
tions in I Zw 18, we have selected three regions (A, B, 
and C) and two regions (E and D) in the main and sec- 
ondary bodies, respectively. In the main body, region A 
is the most crowded one and corresponds to what has 
been referred to as the NW component in previo us pa- 
pers (see e.g. Fig. 1 in IPapaderos fc Ostl in 2012), while 
region C is the least crowded one and is the most suit- 
able to study the old stellar population. In the secondary 
body, region E follows the comma-shaped distribution of 
the brightest blue stars, while region D comprises the 
redder less crowded "halo" . 

The first important result of this paper is an estimate of 
the stellar mass locked-up in old (age>l Gyr) stars. Our 
study indicates a mass of « 1O 6 A7 in old stars in region 
C, although with large uncertainties. This is about half 
of the total stellar mass in region C. On the other hand, 
the severe crowding in regions A and B prevents to reach 
RGB stars and to derive the SFH at ages older than 1 
Gyr. From this result we can derive a lower limit for the 
total mass in old stars in the whole main body. Assuming 
in region A and B the same surface density of old stars 
than in region C, and using the areas reported in Table[TJ 
we obtain a mass of « 2 x 1O 6 M . Since stars tend to be 
more concentrated toward the center, it is likely that the 
true mass in old stars in the main body is significantly 
higher than this. As for the secondary body, the amount 
of old stars is less constrained and we don't attempt to 
quantify it. 

The second important result of our study concerns the 
very recent SF in I Zw 18. The main body has been form- 
ing stars very actively in recent epochs, with an average 
SFR over the last 10 Myr as high as « XM^j-yr in the 
NW component (region A). This corresponds to a specific 
rate of ~ 2x 1Q~ 5 Mqut -1 pc~ 2 , which is just a few factors 
lower than the rate derived in the strong star forming re- 
gion 30 Doradus in the LMC. We recall however that this 
result relies on the assumption that all the brightest stars 
in the CMD are in fact individual objects rather than 
blends of two or more stars or unresolved star clusters, al- 
though we have demonstrated that the colors of the red- 
dest (y — 7 > 1) bri ght objects are incompat i ble w ith the 
cluster hypothesis ( Contreras Ramos et al.l [201 ID . In- 
deed, our value can be compa red with the rate inf erred 
from the Ha luminosity. From [Cannon et "ail ([2002D and 
assuming a distance of 18 Mpc, we derive a Ha luminos- 
ity of ~ 5.9 x 10 39 erg/s for the NW component. Notice 
that this only includes the ionized gas around the peak 
of the continuum emission (wregion A) but not the ion- 
ized gas extendin g outside the ste l lar co mponent. Using 
the relation from iKennicutt et al.l (|1994D for a Salpeter's 
IMF in the 0.1-100 M mass interval, we obtain a SFR of 
^0.05 M /yr, significantly lower than what was derived 
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from the synthetic CMDs. We notice that the occurrence 
of star formation at very recent epochs in the main body 
is in agre ement with the presence of ultra long pe riod 
variables (jFiorentino et al.ll2010fc Marconi et al. 20l3) . A 
potential rate of ~ IM^j/yr in region A implies a stellar 
mass of « 10 7 M Q . This high conce ntration of mass i s 
interesting in view of the results bv iLelli et al.1 |2012), 
who analyzed archival VLA data and found that the HI 
associated to the starburst region forms a compact fast- 
rotating disk. The rotation curve is flat with a steep rise 
in the inner parts, indicating the presence of a strong 
central concentration of mass. 

In the main body, the average SFR over the last 10 Myr 
decreases from the most crowded region A to the least 
crowded region C down to w 2 x 10~ 3 M Q /yr. However, 
this is still comparable to the average rate inferred in 
region C over the last ~ 1 Gyr. On the other hand, the 
secondary body was much less active than the main body 
during the last ~10 Myr, in agreement with the absence 
of nebul ar emission everywhere but in its central region 
(see also IPapaderos fe Ostlinl 120121) . The peak of SF is 
found to occur ~15 Myr ago in the more crowded region 
E and ~ 250 Myr ago in the less crowded region D. 

7. SUMMARY 

• We confirm that I Zw 18 has started forming stars 
earlier than ~1 Gyr ago, and possibly at epochs 
as old as a Hubble time. Thus it is not a truly 
young galaxy at its first bursts of s tar formation, as 
argued in previous studies (e.g.. | Izotov fc Thuanl 
[200irjamet et al.ll20Tot IPapaderos fc QstlirJl2012h . 

• In the periphery of I Zw 18' s main body, where 
crowding is low enough to potentially detect the 
old population, we estimate a mass of « 10 Mq 
in stars older than ~1 Gyr, accounting for about 
half of the total stellar mass in this region. On 
the other hand, crowding is too severe to allow the 
detection of such a population in the more central 
regions. Assuming there the same surface density 
of old stars than in the periphery, we obtain a lower 
limit of s 2 x 10 6 Mq for the total mass in old 
stars in I Zw 18' s main body, The presence of 
a significant amount of old stars in the secondary 
body is more uncertain. 

• The main body has been forming stars at a 
very high rate in recent epochs. In the most 
crowded NW region (region A), the average 
rate over the last ~10 Myr was as high as » 
lMts/yr, corresponding to a specific rate of w 
2 x 10~ 5 MqU^ 1 pc~ 2 . However, this result relies 
on the assumption that all the brightest stars in 
the CMD are in fact individual objects rather than 
blends of two or more stars or unresolved star clus- 
ters. 

• While in the main body the peak of activity oc- 
curred during the last ~10 Myr, the secondary 
body was much less active at these epochs, in agree- 
ment with the absence of significant nebular emis- 
sion. In the secondary body the peak of activity is 
found to occur ~15 Myr ago in the more crowded 
central region and ~ 250 Myr ago in the periphery. 



• The high current SFR in I Zw 18 explains why 
this galaxy is so blue and has a high ionized gas 
content, resembling primeval galaxies in the early 
Universe. Detailed chemical evolution models 
are required to quantitatively check whether the 
SFH from the synthetic CMDs can explain the 
low measure d element abundan c es, both in the 
ioniz ed (e.g., iGarnett et al.lfl997t llzotov fc Thuanl 
1998) and in the n eutral interstellar medium 
( Aloi si et al.1 120031: iLecavelier des Etangs et al.l 
120041 : Uames et al.ll2013l ). or if galactic winds with 
loss of metals are nee de d, as or iginally s u ggeste d 
bv iMatteucci fc Tosil (fl985h ; iPilvuginl (fl993l) : 
iMarconi et al.l (|1994f ). 
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